Effects of patch size and number within a simple model of patchy colloids 
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We report on a computer simulation and integral equation study of a simple model of patchy 
spheres, each of whose surfaces is decorated with two opposite attractive caps, as a function of the 
fraction x of covered attractive surface. The simple model explored — the two-patch Kern-Frenkel 
model — interpolates between a square-well and a hard-sphere potential on changing the coverage 
X- We show that integral equation theory provides quantitative predictions in the entire explored 
region of temperatures and densities from the square- well limit x = 1-0 down to x ~ 0.6. For smaller 
X, good numerical convergence of the equations is achieved only at temperatures larger than the gas- 
liquid critical point, where however integral equation theory provides a complete description of the 
angular dependence. These results are contrasted with those for the one-patch case. We investigate 
the remaining region of coverage via numerical simulation and show how the gas-liquid critical point 
moves to smaller densities and temperatures on decreasing x- Below x ~ 0.3, crystallization prevents 
the possibility of observing the evolution of the line of critical points, providing the angular analog 
of the disappearance of the liquid as an equilibrium phase on decreasing the range for spherical 
potentials. Finally, we show that the stable ordered phase evolves on decreasing x from a three- 
dimensional crystal of interconnected planes to a two-dimensional independent-planes structure to 
a one-dimensional fluid of chains when the one-bond-per-patch limit is eventually reached. 
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I. INTRODUCTION 

Spherically symmetric potentials have become a well-established paradigm of colloidal science in past decades. [l| 
This is because, at a sufficiently coarse-grained level, colloidal surface composition can be regarded as uniform with a 
good degree of confidence, so that relevant interactions depend only on relative distances among the particles. Recent 
advances in chemical particle synthesis have however challenged this view by emphasizing the fundamental role of 
surface colloidal heterogeneities and their detailed chemical compositions. This is particularly true for an important 
subclass of colloidal systems, namely proteins, where the presence of anisotropic interactions cannot be neglected, 
even at the minimal level. [JQ Directional interactions introduce novel properties in such systems. These properties 
depend both on the number of contacts {i.e., the valency) and the amplitude of these interactions {i.e., flexibility of 
the bonds), a notable example of this class being hydrogen-bond interactions, ubiquitous in biological, chemical, and 
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physical processes. [1, Q 

As a reasonable compromise between the high complexity of interactions governing the above systems and the 
necessary simplicity required for a minimal model, patchy-sphere models stand out for their remarkable success in 
this rapidly evolving field. [8l4l]| See Ref. [l2| for a recent review on the subject. 

Within this class of models, interactions are spread over a limited part of the surface, either concentrated over a 
number of pointlike spots or distributed over one or more extended regions. HBl While the former have 

the considerable advantage of a simple theoretical scheme [T6j which allows a first semi-quantitative description, the 
latter can easily account for both the effect of the number of contacts and their amplitude, unlike "spotty" interactions 
which are always limited by the one-bond-per-site constraint. 

In this paper we consider a particular model due to Kern and Frenkel |15| of this patchy-spheres class wherein short- 
range attractive interactions — of the square-well (SW) form — are distributed over circular patches on otherwise 
hard spheres (HS). Interactions between particles (spheres) are then attractive in the SW-SW interfacial geometry 
or purely hard-sphere repulsive under the HS-SW or HS-HS interfacial geometries, and can sustain more than one 
bond — in fact, as many as the geometry allows — even in the case of a single patch assigned to each sphere. A 
number of real systems ranging from surfactants to globular proteins can be described with simplified interactions 
of these particular forms, with well-defined solvophilic and solvophobic regions, and despite their simplicity patchy 
hard spheres have already shown a remarkable richness of theoretical predictions. [ij, [l5|, [T7l - [l9| Notwithstan ding 
the discontinous nature of the angular interactions, highly simplified integral equation approaches are possible, [l7[ 
but only very recently has a complete well-defined scheme, within the framework of the reference hypernetted-chain 
(RHNC) integral equation, been proposed and solved for patchy spheres. [1^ This integral equation belongs to a class 
of approximate closures which have been extensively exploited in the field of molecular associating fluids. [2l| Its 
main advantage over other available approximations (other than its less-accurate parent HNC closure) lies in the fact 
that it relies on a single approximation, for the bridge function appearing in the exact relation between pair potential 
and pair distribution function (?(12), [2ll [22j to directly yield structural and thermodynam ic p roperties that include 
the Helmholtz free energy and the chemical potential with no further approximations. [23l . |24| In addition, it can be 
made to display enhanced consistency among different thermodynamic routes, psj This is an important point when 
analyzing fluid-fluid phase diagrams such as we propose to do here. We thus build upon our previous work with the 
one-patch potential [20| to study the two-patch case and its relationship with its one-patch counterpart. In addition 
to RHNC integral equation results, we provide dedicated Monte Carlo simulations which can assess the performance 
of RHNC. We find that RHNC provides a robust representation of both structural and thermophysical properties of 
the two-patch Kern-Frenkel model for a wide range of coverage x (the ratio between attractive and total hard-sphere 
surface), extending from an isotropic SW to a bare HS potential. The competition arising between phase separation 
and polymerisation is discussed in terms of the angular dependence of the pair correlation function and the structure 
factor. Finally, a comparison between the one-patch and two-patch phase diagrams shows a strong impact on the 
different morphology and stable structures obtained in the two cases. 

We also report numerical simulation results of the model in the region where the RHNC integral equations do 
not numerically converge, to explore the low temperature, small x limit. We find that for x < 0.3 it becomes 
impossible to investigate the low-temperature disordered phases, since the system quickly transforms into an ordered 
structure, which itself depends on the coverage value. Indeed, on decreasing x one progressively enters the region 
where the maximum number of contacts per patch evolves from four to two and eventually reaches the one-bond- 
per-patch condition. When three or four bonds per patch are possible, the observed ordered structure is a crystal of 
interconnected planes, while when only two contacts are possible, particles order themselves into a set of disconnected 
planes. 

The patchy interaction model examined here can be regarded as a prototype of a special colloidal architecture where 
there exist competitive interactions on the colloidal surface that drive, by free energy minimization, the different 
colloidal particles through a spontaneous self-assembly process into complex superstructures whose final target can 
be experimentally probed and properly tuned. [1^ The possibility, discussed in the present study, of identifying the 
position of the gas- liquid coexisting lines and its relative interplay with different structures, opens up fascinating 
scenarios in material science, on the possibility of novel material design exploiting a bottom-up process not requiring 
human intervention. 



II. THE TWO-PATCH KERN-FRENKEL MODEL 

As a paradigmatic model for highly anisotropic interactions, we take the Kern-Frenkcl [ij two-patch model where 
two attractive patches are symmetrically arranged as polar caps on a hard sphere of diameter a. Each patch can 
be reckoned as the intersection of a spherical shell with a cone of semi-amplitude and vertex at the center of the 
sphere. Consider spheres 1 and 2 and let f 12 be the direction joining the two sphere centers, pointing from sphere 1 
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to sphere 2 (see Fig. [T]) . The orientation of sphere i is defined by a unit vector = n^*'' passing outward through the 
center of one of its patches, to be arbitrarily designated as the "top" (i) patch. The patch on the opposite, "bottom" 
(6) pole is then identified with the outward normal hf'-' = hi. 

Two spheres attract via a square- well potential of range Xa and depth e if any combination of the two patches on 
each sphere are within a solid angle defined by 9q and otherwise repel each other as hard spheres. The pair potential 
then reads llSl 



where 



$(12) = 0(ri2)*(ni,n2,fi2), (1) 



oo, < r < a 

0(r) = •( -e, (T <r < Xa (2) 
0, Act < r 



and 



* (ni,n2,fi2) = ( ^' "i*"^ ■ - ^^^'^ ^"2^'^ ■ ^ cos6'o 

' ' 10, otherwise 



(3) 



localized interactions having the single-bond, single-site limitation. [10|, Il3|, |27 
examined potentials of the Kern-Frenkel form using numerical simulations. 



where pi,p2 — ^ or b indicates which patch, top or bottom, is involved on each sphere. The unit vectors ni{iOi) are 
defined by the spherical angles tui = (Oi, tpi) in an arbitrarily oriented coordinate frame and fi2(fi) is identified by the 
spherical angle 57 in the same frame. Reduced units, temperature T* = ksT/e and density p* = pcr^, will be used 
throughout. 

This model was introduced by Kern and Frenkel, [ist patterned after a similar model studied by Chapman et at, 
[l^ as a minimal model where both the distributions and the sizes of attractive surface regions on particles can be 
tuned. In this sense, the model constitutes a useful paradigm lying between spherically symmetric models which 
do not capture the specificity of surface groups, not even at the sim ples t possible level, and models with highly 

y Several previous studies have already 
I, llSj corresponding-state arguments, 

18| highly simplified integral equation theories, and pertubation theories. [19j More recently, [2^ the single 
patch Kern-Frenkel potential was studied using a more sophisticated integral equation approach based on the RHNC 
approximation coupled with rather precise and extensive Monte Carlo simulations. In the present paper, we extend this 
last study to the two-patch Kern-Frenkel potential and provide new methodologies specific for the angular distribution 
analysis. 

We define the coverage x as the fraction of the total sphere surface covered by attractive patches. Thus x = 1 
corresponds to a fully symmetric square-well potential while x = corresponds to a hard-sphere interaction and the 
model smoothly interpolates between these two extremes in the intermediate cases < x < 1. The two-patch potential 
is expected to present qualitative as well as quantitative differences with respect to its one-patch counterpart. One 
interesting question, for instance, concerns the subtle interplay between distribution and size of the attractive patches 
on the fluid-fluid phase separation diagram. It is now well estalDlished 0, [H, [13, H^l that as coverage decreases the fluid- 
fluid coexistence line progressively diminishes in width and height. Indeed, this feature can be exploited to suppress 
phase separation altogether to enhance the possibility of studying glassy behavior [l^, [3 and cannot be accounted 
for with a simple temperature and density rescaling, jTzj although corresponding-state type of arguments can be 
proposed. [Tsj On the other hand, the above mechanism can signiflcantly depend on how the same reduced attractive 
region is distributed on the surface of particles. In Ref. [13, for instance, it was suggested that lines of decreasing 
critical temperature as a function of decreasing coverage, for the one-patch and two-patch Kern-Frenkel models with 
very short-range interactions, could cross each other at a specific coverage: for low coverages, critical temperatures for 
the one-patch model lie above the two-patch counterpart whereas the opposite is true for larger coverages. This would 
have far-reaching consequences on the phase diagram, as phase separation would occur at higher or lower temperatures 
for fixed coverage, depending on the specific allotment of the coverage. Another interesting issue regards micellization 
phenomena, present in the one-patch version of the model, (28j which is expected to be replaced by polymerization 
(or chaining) in the two-patch version. [2^ 



III. INTEGRAL EQUATION WITH RHNC CLOSURE AND MONTE CARLO SIMULATIONS 

The Ornstein-Zernike (OZ) equation defines the direct correlation function c(12) in terms of the pair correlation 
function h(12) = g{12) — 1; it is convenient for computation to write it using the indirect correlation function 



4 



7(12) = h{12) - c(12) instead of h{12). Wc have then 



7(12) ^ £ I dr3dw3[7(13) + c(13)]c(32). 



A second, or "closure," equation coupUng 7(12) and c(12) is needed. The general form for this is [22 

c(12) = cxp[-/3$(12) +7(12) + B(12)] - 1 -7(12), 



(4) 



(5) 



where f3 — (kBT)^^ and a third pair function, the so-called "bridge" function i?(12), has also been introduced. 
While known in a formal sense as a power scries in density, [2^ B{12) cannot in fact be evaluated exactly and at 
this point an approximation is unavoidable. The RHNC approximation replaces the unknown i?(12) with a known 
version _Bo(12) from some "reference" system. In practice, only the hard-sphere model is today well-enough known 
to play the role of reference system. Here we will use the Verlet-Weis-Henderson-Grundke parametrization [sO, IMj 
for i?o(12) = B}is{ri2; ao), where ctq is the reference hard-sphere diameter. Some computational details of the RHNC 
integral equation approach can be found in Ref. [l^ (see expecially Appendix A) , so only the most relevant equations 
will be repeated here. 

Solution of the Ornstcin-Zcrnikc integral equation for molecular fluids [2l| seemingly requires expansions in spherical 
harmonics of the angular dependence of all pair functions, a need that would be very problematic in the case of the 
discontinuous angular dependence in the present $(12). In fact, the integral equation algorithm allows $(12) to 
remain unexpanded. [l^l There is a potential problem however in evaluating the Gauss-Legendre quadratures used in 
the numerical solution, in that of the angles 61,62, .■■ ,0n used for an nth-order quadrature, none is likely to coincide 
with the angle defining the semi-amplitude of a patch. Thus the algorithm will not "know" the correct patch size. 
This problem is ameliorated in the following ad hoc fashion. 

From the interaction $ (12) of Eq. the total coverage x can be computed in terms of 60 as 



1 



X 



{Any 



dLUldLU2 



- COS 6*0) 
COS 6*0) 



(6) 



Q (cos 6*1 — cos(?o) (^ cos 02 

-1-0 (cos 01 — COS 6*0) (cos 6*2 - 

+Q {— COS 61 — COS 6o)Q {— cos 62 — cos 60 

-t-0 (— COS01 — cos^o) (cos6'2 — COS^o) , 

where QJx) is the Heaviside step function, equal to 1 if 2; > and if .t < 0. The integrals can be readily evaluated 
to give [1^ 



= 2 sin 



2 '^o 



(7) 



This quantity can also be numerically evaluated by Gauss-Legendre quadrature using the n roots 6j of the Legendre 
polynomial Pn (cos 6) and the computed result compared with the exact value ([7]) . We may then vary n so as to find 
that number n (typically kept between 30 and 40) that minimizes the known error in computing x- All Gaussian 
quadratures for that x value arc then evaluated with the same number n of points, thus ensuring that minimal error 
arises from the selected angular grid. 

In an axial r frame [2l| with fi2 = z, the internal energy per particle in units ksT is obtained from 



N 



dr {g{r,uJi,uJ2) * (wi, a;2))a;ia 



(8) 



where (. . = (l/47r) J duj . . . denotes an average over spherical angle w and where we have written out ^(12) — 
g{r,LjJi,uj2). Similarly, the pressure P is computed from the compressibility factor 



pP_ 
P 



+ (y(fT,Wi,c.2)e^^*("^'"=)\ -A3/y(Aa,wi,a;2) 

O I \ / UJ\UJ2 ^ 



(9) 



where the cavity function 2/(12) = g{12)e^^'^^'^^ has been introduced. The angular integrations in these expressions are 
evaluated with Gauss-Legendre and Gauss- Chebyshev quadratures. Finally, the dimensionless free energy per particle 
/3F/N and chemical potential (3fj, can also be directly computed from the pair functions produced by the RHNC 
equation; the overall calculation is optimized by choosing the reference hard sphere diameter co so as to minimize the 
free energy functional. [2^ We solve the RHNC equations numerically on r and k grids of A^^ = 2048 points, with 
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intervals Ar = O.Olcr and Afc = 7r/(iVrAr), using a standard Picard iteration method. [23| The square-weh width is 



set at A = 1.5 as a reasonable value dictated by the availability of isotropic square well results. [32| Further details of 
these and other computations can be found in Ref. [IS- 

For an assessment of the performance of the RHNC integral equation, we also perform NVT, grand canonical, and 



Gibbs ensemble Monte Carlo (MC) simulations [33[ following the path set in the one-patch case. [20| Standard NVT 
MC simulations of a system of 1000 particles are used to compute structural information (pair correlation functions 
and structure factors) for comparison with integral equations results, whereas grand canonical and Gibbs ensemble 
MC (GEMC) are used to locate critical parameters and coexisting phases. The exact locations of the critical points 
(points connected by the thick dashed green line in Fig. [2]) have been obtained from the MC data assuming the Ising 
universality class and properly matching the density fluctuations with the known fluctuations of the magnetization 
close to the Ising critical point. [U For GEMC, we use a system of 1200 particles, which partition themselves into 
two boxes whose total volume is 4300cr'^, corresponding to an average density of p* = 0.27. At the lowest temperature 
considered, this corresponds to roughly 1050 particles in the liquid box and 150 particles in the gas box (of side sa IZa). 
On average, the code attempts one volume change every five particle-swap moves and 500 displacement moves. Each 
displacement move is composed of a simultaneous random translation of the particle center (uniformly distributed 
between ±0.05(t) and a rotation (with an angle uniformly distributed between ±0.1 radians) around a random axis. 
We have studied systems of size L = 7 up to L = 10 to estimate the size dependence of the critical point, with 
an average of one insertion/deletion step every 500 displacement steps in the case of grand canonical Monte Carlo 
(GCMC). We have also performed a set of GCMC simulations for different choices of T and to evaluate p*{fi,T). 
See Ref. [13 and additional references therein for details. 

IV. NUMERICAL RESULTS 
A. Coexistence line 

Locating coexistence lines is not an easy task within integral equation theory, given the fact that virtually all 
integral equations are unable to access the critical region with reliable precision due to significant thermodynamic 
inconsistencies among various possible routes to thermodynamics, a consequence of the approximation buried in the 
closure Eq. ([S]). The RHNC closure is no exception to this rule, but has the strong advantage of relying on a single 
approximation expressed by the choice of the reference bridge function i?o(12), at odds with other available closures 
which require additional approximations in constructing various thermodynamic quantities such as the chemical 
potential. Here we follow the protocol outlined in Ref. [32| for the isotropic square- well potential and Ref. [13 for the 
one-patch Kern-Frenkel potential, where both the well-known pseudo-solutions shortcoming [ssj and the numerical 
drawbacks [36| can be conveniently accounted for. 

Figure [5] depicts the location of the fluid-fluid coexistence line for the two-patch case upon varying the coverage x- 
The limiting case x = 1 corresponds to the square-well potential. Both MC (points) and RHNC (thick solid lines) 
results are shown. As previously noted, RHNC is not able to approach the critical point close enough to provide a 
direct estimate of its location. However, since it provides a quite good description of the low-temperature part of the 
coexistence line, we have tried to use these data to approximately locate the critical point. 

Visual inspection of the RHNC coexistence points reveals, in the cases where it is possible to go closer to the 
critical region, an unphysical change of curvature of the coexistence line moving from low to high temperature. For 
this reason, for each coverage we selected only data clearly consistent with a rectilinear diameter law. Then we fitted 
those data with the following function (corresponding to the first correction to the scaling (STj): 

PI- Pg = a{T,-Tf{l + b{T,-T)^), (10) 

where the values of the exponents /3 = 0.325 and A — 0.54 are appropriate for the 3D Ising model universality 
class. [3^ [S^l Once Tc and the amplitudes a, b have been determined, the critical density can be obtained from 
the rectilinear diameter best fit. The numerical results of such a procedure arc compared with MC estimates of the 
critical points in Table [H It is evident that even though in general the validity of Eq. ((TU)) is deemed to be limited to 



a smaller neighborhood of the critical point, [38| in the present case it provides an acceptable procedure for a quick 
first estimate of the critical point location. 

As the coverage decreases, the coexistence line shrinks and moves to lower temperature and density, as expected 
from an overall-decreasing attractive interaction. This trend can be tracked rather precisely by MC simulations down 
to remarkably low coverages (x = 0.3) and RHNC correctly reproduces this evolution down to x = 0.6 coverage. 
Below this value, more powerful algorithms are required to achieve good numerical convergence. 

A few remarks are here in order. The coexistence curves shown in Fig. [5] arc consistent with previous analogous 
results reported in Ref. [l^ but extend the range of temperatures and, more importantly, the range of coverages (x 
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values). This allows a quantitative measure of the significant deviation from the simple mean- field-like results which 
can be obtained from the simple scaling (not shown) T* T* / x a,s suggested by the second-virial coefficient B2{T*) 
for this model, [lH 

= 1 - x^A^ - 1) (^^^ - 1) - (11) 

B2 being the hard-sphere result. This is also consistent with the breakdown of the above simple scaling at the level 
of the third virial coefficient, derived in Ref. [l3 for the companion patchy sticky-hard-sphcre model. As we shall see 
later on, the dependence of the critical temperature and density on both the coverage and the number of patches is 
one of the main results of the present work. As a final point, we note that all curves in Fig. [5] collapse into a single 
master curve upon scaling T* — > T* /T*. in agreement with Ref. ITsl 



B. Low-coverage results 



Below X — 0.3, it becomes impossible to properly estimate the location of the critical point or the density of the 
coexisting gas and liquid phases. Indeed, the gas-liquid separation becomes pre-empted by crystallization into a 
structure that depends on the value of x- Hence, the liquid phase, as an equilibrium phase, ceases to exist for small 
X- This is strongly reminiscent of the disappearence of the liquid phase using spherical potentials when the range of 
the interaction becomes smaller than about 10% of the particle diameter, [40l - |42| thus providing the angular analog 
of the same phenomenon. Interestingly enough, the crystal structure which is spontaneously observed during the 
simulation depends on the value of x, since the x value controls the maximum number of bonds per patch. In the 
range of x values such that each patch can be involved in four bonds, the observed ordered structure is made by planes 
exposing the SW parts to their surfaces (see Fig. [Sl[d)). Particles in the plane are located on a square lattice and 
adjacent planes are shifted in each direction by a half lattice constant, resulting in a reduced energy of -4 per particle 
(i.e., eight bonded neighbors). On decreasing x below x ~ 0.118, the region where only three bonds per particle are 
possible ([•\/3(l-|- A/cr)]~^ < sin^o < W^{^ + is entered and the crystal structure is made by interconnected 

planes of particles arranged in a triangular lattice (see Fig. [3Jc)). For [a/4(1 -I- A/cr)]~^ < sin^o < [a/3(1 + 
only two bonds per patch are possible and the system organizes into independent planes (see Fig. [Sjb)), this time 
turning a HS surface to their neighboring planes. Particles in the plane are now arranged on a triangular lattice and 
each patch is able to bind only to two different neighbors located in the same plane, resulting in a reduced energy 
per particle of —2. When sin^o becomes smaller than the value [2(1 -I- A/ct)]~^ (corresponding to x ~ 0.0572), the 
one-bond-per-patch condition is reached and the system can form only isolated chains (see Fig. [3l^a)). In this limit, 
the system is expected to behave as the two single-bond-per-patch model. j29| 



C. Structural information 



We turn our attention next to structural information, where the advantages of a reliable integral equation approach 
become evident. One has to keep in mind that Gibbs ensemble and GCMC simulations are particularly painstaking, 
due to the combined effect of the required low temperatures and the aggregation properties of the fluid (as detailed 
below), so that many of the state points examined here require several weeks of computer time. On the other hand, 
the RHNC integral equation, while rather demanding from an algorithmic point of view (see e.g. Appendix A of 
Ref. [20I ) is a rapidly convergent scheme yielding solutions on the order of minutes, depending on the temperatures 
considered. A more profound advantage stems from the fact that, within the approximation defined by the RHNC 
closure, all possible pair structural information is in fact exactly available, unlike MC calculations where, though 
available in principle, their statistics would be so limited as to make such calculations impractical. Thus, only the 
pair correlation function g(12) averaged over angle Yi2{VL) is computed. By symmetry, the resulting pair function in 
this context depends only on r = ri2 and cos0i2 = fii • n2 and will be denoted here as g(r, cos (^12); see Appendix B 
in Ref. ^ for details. 

Consider then the unaveraged pair correlation function .g(12) = (7(r, wi, 012) from RHNC in an axial frame with 
fi2 = z. Two noteworthy configurations occur when (a) all four patches lie along the same line (we denote this as 
the parallel (||) configuration, with fii • 112 = ±1, the actual labeling of each patch being unimportant) and when (b) 
patches on sphere 2 lie on an axis perpendicular to those of sphere 1 (in this case ni • n2 = 0, which we denote as the 
crossed (X) configuration). 

Figure m reports the results for the case A = 1.5, p* = 0.7, and T* = 1.0, which has been selected so that the fluid 
is above the coexistence line for all considered coverages, with conflgurations || and X in the top and bottom panels 
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respectively. Values of coverages range froni a full square-well potential (x = 1) to a hard-sphere potential (x = 0). 

As coverage decreases, the contact value r = cr"*" of the || configuration has the unusual behavior of first a slight 
increase from x = 1 to x = 0.5, followed by a more marked increase starting at x = 0.4 up to the very small coverage 
X = 0.1 limit which eventually backtracks to roughly the same value as at x = 0-4 in the hard-sphere limit. At 
the opposite side of the well, r = (Act)^, an even more erratic behavior is observed, with an increase in the range 
0.7 < X < 1; then a decrease for 0.3 < x < 0.6, a new increase down to x = 0.1, and a final sudden decrease to 
the hard sphere value x = 0. A somewhat similar feature occurs in the X configuration where within the entire well 
region cr"*" < r < (A(t)~ one observes a sudden decrease from x = 1 to x = 0-9 ^-^d a more gradual increase until 
reaching the highest value for the HS case. It is worth noting that in the X configuration there is no discontinuous 
jump at r = Act for any value x < 1- The reason for this has already been addressed in Ref. [13 for the one-patch 
case. Outside the first shell, there is a very weak dependence on the coverage, with a slight shift in the location of 
the second peak from a value of r sa 2.25ct at the SW x = 1 to a value of r « 1.8ct for lower x- 

We compare RHNC integral equation results with MC simulations in Fig. [S] As noted above, only the averaged 
pair function g(r, cos0i2) can be compared and this is done in the figure for different values of the coverage x at the 
same state point (T* = 1.0, p* = 0.7) and for the same || and X configurations considered earlier. The good overall 
performance of RHNC in representing MC results is apparent as both contact values at the well edges and the jump 
discontinuities are very well reproduced. It is instructive to contrast these results with those of Fig. |4l as many of 
the abrupt changes appearing in the actual pair correlation function are smoothed out by the orientational average 
carried out here. For instance, the characteristic jump at r = Act of the 1 1 configuration (top panel) progressively 
decreases as coverage is reduced and disappears in the hard-sphere limit. Conversely, the jump is present also in the 
X configuration (bottom panel) unlike the corresponding case of the full 5(12). In addition, the strong increase of the 
II configuration for low coverages is not present in this figure; as remarked earlier, this level of detail in (7(12) is one 
of the main advantages of an integral equation approach. 

An additional useful quantity to consider, in view of its direct experimental access through scattering experiments, 
[2^ is the structure factor, which will be denoted S'ooo(fc) within our theoretical framework. [13, HH, [ll] This is also 
strongly related, via Hankel transforms, to the radial distribution function .9000 (j'), which is g(12) averaged over all 
orientations wi and uj2 of the patches and of the relative angular position fi. Note that (7ooo('') is also the simplest 
rotational invariant (see Ref. [2l| and Appendix |B]) . 

The structure factor and the radial distribution function are reported in Fig. [B] for two representative values of 
coverage, x = 0-8 and x = 0.2, corresponding to almost fully attractive and almost fully repulsive limits. These 
values have been selected at the same state point previously considered (T* = 1.0, p* = 0.7) as having a very different 
behavior within the first shell ct < r < Act. This high-density result is also contrasted with a low-density state point 
p* = 0.1 at the same temperature, a value which, in the temperature-density plane, lies symmetrically with respect 
to the coexistence curves in the single fluid phase for all coverages (see Fig. 

A few features are worth noting. For density p* = 0.7 there is a significant coverage dependence, where the contact 
value 3000(0"^) for x = 0.2 coverage is larger than that for the corresponding x = 0.8 case and, conversely, the jump 
present at the other extreme Act~ is much smaller in the former than in the latter case. A similar feature also occurs 
for the low-density state point p* =0.1. This results from an angular average of the results given in Fig. 21 Likewise, 
there is a marked difference in the behavior of the structure factor 6*000 (fc) for the high-density case p* = 0.7, both 
in the height of the first peak (related to the 5000 (f^) value) and of the secondary peaks (related to the behaviors of 
5000 ('') in the ct < r < Act region and of the gooo(-^f ^) discontinuity). Similarly, in the low-density branch p* = 0.1, 
the large 6*000(0) value for the x = 0.8 coverage case is signaling the approach to a spinodal instability which is clearly 
not present in the corresponding x = 0.2 coverage. 

One natural interpretation of the above results is the progressive rearrangement of the distribution within the first 
shell upon varying both the coverage and the density. To support this view, we consider the angular distribution 
within the first shell in the next subsection. 



D. Angular distribution 

The nonmonotonic dependence of g{l2) in terms of the distance r/a for decreasing coverage x, as illustrated in Fig. 
131 is rather intriguing and requires an explanation. A similar, albeit different, feature occurs even in the one-patch 
case, as shown in Ref. [l^- We have tackled this in two ways, illustrated in the following. 

All previous representations of g{12) have been depicted in the molecular axial frame, where rii • fi2 = 1, so that 
patches on sphere 1 are parallel to the vector ri2 joining sphere 1 with sphere 2. This is clearly preventing an 
understanding of the angular distribution of the patches around a given sphere 1, that is, as a function of fi2(r2). 

This is however a needless restriction, as one can start from the expression for 5(12) in a general (laboratory) frame 
and study the dependence on the angle for fixed patch directions fii and n2. In Figs. |31|Sland|Hl we notice that the 
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main dependence on coverage x stems from the region within the well, t < r < Act; it sufScies therefore to investigate 
the average f2 dependence by integrating over the radial variable r within this region. We further note that there is 
azimuthal symmetry with respect to the variable, so that we can focus on the dependence. The details of the 
analysis are reported in Appendix El where it is shown that the relevant quantity is g{0,02), which is a function of 
the angle 6 (polar dependence of fi2(r2)) and of the polar angle 02 of the patches on particle 2, given that the patches 
on particle 1 lie along the z axis. We report comparative calculations for both low {p* = 0.1) and high {p* = 0.7) 
density at identical temperature T* = 1.0 at two representative coverages, x = 0-8, representing a case with almost 
all attraction, and x = 0.2, as representative of an almost hard-sphere case. These are the same conditions considered 
in Fig. ini the results are reported in Fig. [T] Let us consider first the high-density, p* = 0.7, situation as depicted in 
the two left panels for x = 0.8 (a) and x = 0-2 (c). Here, the x = 0-8 case yields a very well-defined pattern with a 
periodically modulated distribution of the patches in symmetrical fashion as indicated by the trimodal distribution as 
a function of the relative positional angle 6, so that 0,7r/2,7r are almost equally represented. (Note that 9 = 0,tt are 
necessarily equivalent due to the up-down symmetry of the two-patch distribution.) The two interstitial minima are a 
consequence of the reduced valency — the corresponding fully symmetrical result under this condition would be a flat 
distribution around the value 1.66 in between the two maxima and minima — so this slightly favours perpendicular 
orientation of the patches along the forward (or backward) direction and parallel orientation along the transversal 
direction. Under low coverage (x — 0.2) conditions, on the other hand, there is clear evidence of a parallel orientation 
of the patches along the forward (or backward) direction, the opposite being true for a perpendicular orientation of 
the patches. This confirms the tendency to filament formations previously alluded to. The situation is even more 
evident at low density, p* =0.1, as shown in the right two panels (b) and (d). 

E. Coefficients of rotational invariants 

Additional insights on the angular correlations of patch distributions can be obtained by considering other coef- 
ficients ft,'i'^'(r) = g'i'^'(r) — Si-^i^ifioo of the rotational invariants ?/j'i'^'(r) as defined in Eqs. (|A1[) and (jA2p : they 
have proven to be of invaluable help in discrimina ting between parallel and antiparallel configurations occurring in 
different models such as dipolar hard spheres |44| | and Heisenberg spin fluids. [1^ Some relevant properties of 
these coefficients are also listed in Appendix [B] where we display explicit expressions for the first few coefficients. 

In the two-patch case, we note that all coefficients with h or I2 odd vanish, so we depict the first nonvanishing 
coefficients /i^^"(r), /i^^^(r), h^'^^{r) ~ h'^^^{r) in Fig. [5] for the same state points as before. Note that the left- most 
two curves, (a) and (c), correspond to density p* = 0.7, temperature T* = 1.0, coverages x = 0-8 (a) and x = 0.2 
(c), and are plotted on the same scale. While qualitative trends are similar, the two cases have significantly different 
behavior. Within a given coverange h?'^^{r) and /i^^^(r) are almost coincident, with positive correlation in the well 
region cr < r < Ar, whereas hP'^'^ir) has decreasing positive correlation in the same region and negative correlation for 
r > Act. Numerical values, on the other hand, differ among each other, with small values for h?^^{r) and /i^^^(r) in 
the high coverage case x = 0.8 and significantly higher values in the low coverage limit x = 0.2. 

Likewise, for the low-density state point p* = 0.1, T* = 1.0, right-hand-side plots (b) and (d) can be unambiguously 
discriminated between high (b) and low (d) coverages. We shall return to this point in the comparison with the one- 
patch results, where the physical meaning of these coefficients will be discussed. 

V. COMPARISON WITH ONE-PATCH RESULTS 
A. Phase diagram 

Figure [S] shows the critical parameters for the case of particles with one and two patches, both reported as a function 
of the total coverage. Here only MC results are reported in view of their precision and reliability. The x = 1 limit 
corresponds to the SW case and coincides for both models. An analogous figure, for the case of adhesive patchy 
spheres (the limit of the present model for vanishing ranges), has been reported in Ref. [TtI within a simplified integral 
equation scheme. 

With respect to the adhesive limit, the range of coverages which can be explored numerically is significantly wider 
(for both one-patch and two-patch cases). In the case of two patches, crystallization pre-empts the possibility of 
exploring the smaller values. In the one-patch case, the process of micelle formation, also observed experimentally, 
[46| suppresses the phase-separation process at small x values. 

The critical parameters decrease on decreasing x for both one-patch and two-patch models. The behavior of Tc 
can be explained on the basis of a progressive reduction of the attractive surface. The decrease in the critical density 
becomes significantly pronounced only for the smallest x values and can be attributed to the lower local density 
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requested for extensive bondin g;. S uch behavior is analogous to the suppression of the critical density observed when 
the particle valence decreases. jlOl . l47j In the x region where it is possible to evaluate the critical parameters, Tc and 
Pc for the two-patch case are always larger than the corresponding one-patch values, a trend which can be tentatively 
rationalized on the basis of the ability to form a larger number of contacts and higher local bonded densities for the 
case in which both poles of the particles can interact attractively with their neighbors. No evidence of a crossing 
between the two geometries is observed. Such a crossing has been predicted by a theoretical approach based on a 
virial expansion up to third order in density and appropriate closures of the direct correlation function. [TtI 

It is worth emphasizing that the above dependence on the number of patches, at a given coverage, provides clear 
evidence of the impossibility of rationalizing the change of the critical line on the basis of a trivial decrease of the 
attractive strength of interactions due to the reduction in coverage, as alluded to in Section HV Al 

For the sake of completness, we also report RHNC integral equation results for the most relevant thermodynamic 
quantities, as a function of the coverage x- These are shown in Table lU for the same high-density state, p* = 0.7 at 
T* = 1.0, considered above for structural information. Here we present the reduced internal energy per particle and 
the reduced excess free energy per particle, U /Ne and pPc^/N, respectively, the reduced chemical potential f3p, the 
compressibility factor 0P/p, and the inverse compressibility /3 (dP/dp)rp. These results may be compared with those 
of Table IV in Ref. |20| listing the same quantities for the one-patch counterpart. (We ignore the tiny difference in 
densities between the two calculated states.) The last two columns give the reference HS diameter cto stemming from 
the variational RHNC scheme (see Ref. [13 for details) and the average coordination number within the wells z, whose 
one-patch counterparts are included in Tables IV and V of Ref. [13, respectively. Note that z here is systematically 
larger than in the one-patch case, in qualitative agreement with the MC results of Fig. ^ 

B. Angular distribution and coefficients of rotational invariants 

Within the RHNC integral approach, the analysis of the angular distribution of patches within the first shell given 
in Scction llV DI revealed that the cylindrical symmetry of a pair of patches (2P case) on each particle was very effective 
in driving the system to morphologically different configurations in the low (20%) and high (80%) coverage limits, 
as illustrated in Fig. [T] It is natural to expect a very different situation in the single-patch case (IP case). This is 
indeed the case as further elaborated below. 

For the single patch with x = 0-2 coverage (Fig. [TOl bottom panels (c) and (d)) parallel patches (6*2 = 0) are more 
likely in the perpendicular direction [0 ^ tt/2 or cos6' w 0), whereas antiparallel patches {O2 — tt) arc more likely in 
the forward direction (cosO « 1). The case of perpendicular patches {62 = ti'/Z) are conversely more or less equally 
distributed along the whole angular region Q < 6 < i:. There is no qualitative difference between the situation of high 
[p* ~ 0.7) and low (p* = 0.1) densities as shown by the contrast between the bottom left (c) and right (d) panels. 
Note that the result significantly contrasts with the corresponding results of the two-patches case depicted in Fig. [T] 
Consider now the opposite situation of very large coverage (x = 0.8) (Fig. [TOl (a)) where there is a single well-defined 
peak for antiparallel orientations [62 = tt) in the backward direction [cosO w — 1). Again, this markedly differs from 
the two-patch case (Fig. [71 top left panel (a)), where there is a triple peak for aligned patches {O2 = 0,7r) in the 
forward {cosO ~ 1), perpendicular (cos 6* w 0), and backward (cos 6' « —1) orientations. This is a dense state point. 
Under diluted conditions, p* = 0.1, we find a qualitatively similar behavior as in the dense case, with antiparallel 
alignment in the forward direction (which cannot be physically distingushed from the backward one). Clearly the 
predominant antiparallel alignment is reflecting the tendency to micellization rather than polymerization which is 
built into the single patch symmetry. 

It is also interesting to contrast the coefficients of rotational invariants for the one-patch case with those obtained 
in the two-patch counterpart in Fig. [71 At variance with the two-patch case, here all coefficients are nonvanishing so 
that we consider the first nonvanishing instances, that is h^^^{r), h^^'^{r), and /i^^°(r), which are particularly useful 
as giving the projections over the important invariants. (43| 

We have evaluated these coefficients for the same state points considered in the two-patch case in Fig. [TT] for 
both dense or diluted conditions and small or large coverages. In contrast with the two-patch case, here the effect of 
coverage appears to be less significant, as can be inferred by inspection of the dense case p* = 0.7 (left panels (a) and 
(c) ). For the ft.^^°(r) case — the projection coefficient along the ferroelectic invariant A(12) in Appendix [B] — we 
find a negative correlation within the well both for x = 0.8 (top left panel (a) ) and x = 0-8 (bottom left panel (c) 
), as expected from the tendency to form antiparallel alignments. Likewise, the projection /i^^^(r) along the dipolar 
invariant D{ni, £12, ri2) is found to be negative and numerically similar to h^^'^(r) at both coverages, again indicating 
negative correlation to dipolar alignment. The only positive correlation is found for the /i^^°(r) component, which 
does not distinguish between up and down symmetry, in qualitative agreement with the two-patch analogue. 

This situation is replicated in the diluted case (right two panels (b) and (d) ) with different numerical values, thus 
indicating that these correlations are signatures of robust orientational trends induced by the the particular one-side 
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symmetry of the single-patch potentiaL 

VI. CONCLUSIONS AND OUTLOOK 

We have performed a detailed study of a fluid whose particles interact via a two-patch Kern-Frenkel potential that 
attributes a negative square-well energy whenever any two patches on the spheres arc within a solid angle associated 
with a predefined coverage and within a given distance given by the well width, and a simple hard-sphere repulsion 
otherwise. This model can be reckoned as a paradigm of a unit system with incompatible elements {e.g., hydrophobic 
and hydrophilic) that can self-assemble into different complex superstructures depending on the parameters of the 
original unit {e.g., coverage). We have exploited state-of-the-art numerical simulations (standard Metropolis, Gibbs 
ensemble, and grand canonical Monte Carlo) coupled with RHNC integral equation theory following the approach 
outlined in previous work on a single patch. |20[ On comparing RHNC integral equation with numerical simulations, 
we find the former to be quantitatively predictive in a large region of coverage, even close to the gas-liquid transition 
critical region, over a range of coverage which is significantly larger than the single patch counterpart. The reason for 
this is attributed to the fact that RHNC uses the approximated hard-sphere bridge function, which retains spherical 
information, as a unique approximation throughout the entire calculation, a feature which works better for the more 
symmetric two-patch case than the highly asymmetric one-patch Kern-Frenkel potential. 

Having assessed the reliability of the RHNC integral-equation approach, we have fully exploited its capabilities 
in providing detailed angular information that is typically inaccessible to MC simulations, as already discussed in 
Ref. [20I . This has been done in two ways. First, by computing the orientational distribution probability of parallel 
and perpendicular alignment of patches within a spherical shell in the region a < r < Act. This methodology is 
able to account for the erratic coverage dependence of the pair correlation function g{12) by clearly discriminating 
between small and large coverages at all densities. The same approach also enlightens the characteristic symmetries 
of the patch distributions when the two-patch case result is contrasted with the one-patch analogue. Second, by 
computing the rotational-invariant coefficients that are the projections of 3(12) over rotational invariants. Again, this 
can discriminate between small and large coverages (at all densities) and single and double patches. 

Our Monte Carlo results extend those originally obtained by Kern and Frenkel [l^ for the two-patch case and 
can be contrasted with those of the corresponding single-patch counterpart [20j and those obtained when the radial 
part of the potential is of the Baxter type. |17| The RHNC calculation presented here, along with the corresponding 
calculation carried out in our previous paper, [20| together constitute the first attempt to apply a well-defined integral 
equation theory to such highly anisotropic potentials having sharp angular modulation. 

An important outcome of our calculations is the clarification of the combined effect that size and distribution of 
the patches have on the gas-liquid coexistence lines and critical parameters. The reduction of the bonding surface 
clearly decreases the critical temperature, an effect which can be related to the decrease in the bonding energy of the 
system. More interestingly, it also shows a suppression of the critical density, which can be interpreted along the same 
lines used in interpreting the valence dependence in patchy colloids. [l3l.l47( Indeed, the maximally bonded structures 
require lower and lower local densities on decreasing x- Interestingly, while in the singlc-bond-per-patch condition 
the evolution of the critical parameters on decreasing valence can be followed down to the limit where clustering 
prevents phase separation, [28| in the model studied here crystallization pre-empts the observation of the liquid-gas 
separation for x < 0.3. Crystallization is here much more effective due to the analogy between the local fully-bonded 
configuration and the crystal structure. By contrast, crystallization is never observed for the one-patch case, where it 
has been shown that the lowest energy configuration is reached instead via the process of formation of large aggregates 
(micelles and vesicles) or via the formation of lamellar phases. (2^ This difference highlights the important coupling 
between the orientational part of the potential and the possibility of forming extended fully-bonded structures. Our 
results indicate that, for a given coverage, in the two-patch case both the critical temperature and density are slightly 
higher then their corresponding one-patch counterparts, thus indicating that an increase of the valence favours the 
gas-liquid transition, in agreement with previous findings. 

A final important consequence of our study concerns the limit of very small coverages that is particularly interesting. 
Indeed, it is possible to tune the structure of the system and control the topology of its ordered arrangement. By doing 
this we have observed a progression from the case where chains are stable (in the one-bond-per-patch limit) to the 
case where independent planes are found, evolving — for slightly larger x values — into an ordered three-dimensional 
crystalline structure. Each of these ordered structures is observed in a restricted range of x values. This possibility of 
fine tuning the morphology by controlling the patterning of the particle surfaces may offer an interesting possibility 
for specific self-assembling structures. 

An additional perspective of our work should be stressed. Several studies (see e.g. Ref. [H and references therein) 
have exploited spherically-symmetric potentials to mimic effective protein-protein interactions, especially in connection 
with protein crystallization. [53 | This is clearly unrealistic for the majority of proteins where the distribution of 
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hydrophobic surface groups is significantly irregular, a feature that can be captured, at the simplest level of description, 
by the model studied here. The specific location of the coexistence lines, such as those considered in the present study, 
have important consequences in the study of pathogenic events for sickle cell anemia [HH and other human diseases. 

M 
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Appendix A: Angular properties of g{12) in a general frame 



The expansion in spherical harmonics Yim{uj) of 5(12) in an arbitrary space frame reads j48l | 

oo h+h 

5(12) ^ ^TT E 9ir;hl2l)i^'''''iuJiL02n), (Al) 

h,l2=0l=\h-l2\ 

where we have introduced the rotational invariants [2l| 

+h +h 

i^'^'-' {LUiu:2^) ^ H C {hhl\mim2mi+m.2)Yi,,nA^i)yh<^^A^2)Yi*^^^^^{n) . (A2) 

mi — — /i 7712 — — ^2 

Note that g{r;lil2l) coincides with ^'^'^'(r) up to a normalization constant (see Appendix IB]) . 

We are free to set the origin of the coordinate frame at the center of particle 1 and choose its orientation so that 
z = ni without loss of generality. We first note that 



2^1 + 1 
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(A3) 



Clearly, the orientation of the x^y axes is then irrelevant, so we may integrate out the angles Lp2 and this leads to 
the average {g{\2)) ^.^^^ where we note that 
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and similarly for {Yi^^ {0:f)) ■ Here the Pm{x) ~ Pi{x) are the usual Lcgendre polynomials. We have then from 
Eqs. (fJT|) and (jM]) that 
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We are interested in the angular behavior within the well, cr < r < Act, so we finally integrate over the radial variable 
r and define 
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drg (r; I1I2I) , 
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The result is then a function of the polar coordinate 6 of fi2 and the polar orientation of the second patch, 62] that 
is, 



(2^1 + 1) (2/2 + 1) (2^ + 1) 
47r 



1/2 



C {hhl^Qm) Pi, (cos 02 ) -Pi (cos 0), 



(A8) 



given that the z axis is aligned with the patches of particle 1. 



Appendix B: Coefficients /i ^ ^ (r) of rotational invariants 

In this Appendix we consider the coefficients ^'^'^'(r) of rotational invariants that have proven to be particularly 
useful in discriminating among different orientational behaviors. In numerical simulations they are defined as follows 
(see e.g. Ref. IH), 
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where the A'i'^'(12) are rotational invariants. Explicit expressions for the first few are [49j 

A°"°(12) = 1, (B2) 

A"° (12) = 3A(12) = 3ni • 62, 

A"2 (12) = Id (12) = ^ [3 (ni • fi2) (n2 • fi2) - ni • fis] , 

A220(12) = ^£;(12) = ^[3(ni.n2)^-l 

Other expressions can be found in Ref. 

We note that the first expression in Eqs. (|B2p yields (7™°(r), which coincides with the radial distribution function 
3ooo('') = (,9(r,wi,aj2))tJiii)2- Here we have h^^^(r) = g°™(r) — 1; in all other cases Ir^^^'^^ir) = g'i'^'(r). Some of the 
coefficients have particularly interesting physical interpretations: the term ft,^^*'(r) is the coefficient of ferroelectric cor- 
relation, the term /i^^^(r) the coefficient of dipolar correlation, the term /i^^°(r) the coefficient of nematic correlation, 
and so on. 

It might be useful to show how these general expressions (typically computed in Monte Carlo simulations) connect 
with the corresponding ones typically evaluated in an integral equation approach. We do this for the representative 
case of g^^^(r\ the others being similar. 

We define ^'^'^'(r) oc g (r; hhl), where the proportionality constant is obtained through a particular prescription to 
be further elaborated below. The projections g{r;lil2l) of g(12) on the rotational invariants '!/''^'^'('^iW2^^) as defined 
by Eq. (jAip are related to the values giii2m{T) in the axial frame by [2l| 

/ 47r V'^ 

gir-Jihl) = [ ^ij^^i ) {Iil2l;mm0)gi^i,„i{r) , (B3) 

where rh = —m. 

Consider as a representative example the quantity A^^^(12) defined in Eqs. (|B2[) and note that in Eq. (jBl[) one 
has 

^ J2 Hr - r,j) D (ij)^ ^ j dujidw2g{r,uJi,U2)D{li). (B4) 

Upon choosing the molecular frame (z = fi2), one finds 

D [12) = 2 cos 01 cos 6*2 — sin 6*1 sin 6*2 cos ((/3i — (/32) (B5) 
47r 

= — [2Yw (wi) Yw (^2) + Yii (wi) yi_i (W2) + yi-i i^i) Yii {LU2)] , 
where the Yimiuj) are spherical harmonics. (2l| 
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The expansion (|A1I) in a molecular frame reduces to 

g{r,uJi,ui) = 47r ^ gij^m {r)Yi^m {uJi)Yi^rii {c^2) ■ 



l\l2m 



Using Eqs. (|B5I) . (|B6I) . and the orthogonality relations for spherical harmonics, [21[ one finds easily that 



1 



(47r)2 



dujiduj2g {r,uJiUJ2) D {12) = [2giio (r) + (r) + gn^i (r)] 



(B6) 



(B7) 



so that, combining with /i^^^(7') given in Eqs. (|Bip and (|B2p and using the symmetry property = .giii(''), 

one finds 



/I'^'W = 3110 W +.9111 (r) 



(B8) 



We can work out the first few projections by using tabulated values for the Clebsch-Gordan coefficients 
C {hhl, mfhO). Using the symmetry properties of the Clebsch-Gordan coefficients and the 5/ii2mi(^)i one finds 
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r;000) = (4^)^/'ffooo (r) , 

r; 110) = - (47r/3)i/2 [gno (r) - 2.gni (r)] , 

r;112) = (8^15)'/' [.9110 W +.9111 (r)], 

r;220) = (MS)'^' [.9220 W - 2,9221 W + 2.9222 W] , 

r;222) ^ - (87r/35)'/' [3220 (r) - 9221 W - 2.9222 W] , 

r;224) = (87r/35)^/^ 

r;011) ^ (4^3)'/' 5010 -.9 (r; 101), 

r;022) = (4^5)'^' 5020 W = g (r; 202) , 

r;121) = -(87r/15)'/' [9i2oW-\/3.9i2i(r)l =-5(r;211). 



(B9) 



4 1 
5220 (f) + ^.9221 (r) + -5222 (r) 



r;123) = (127r/35) 



1/2 



5120 (r) + -^9121 (r) 



-5(r;213). 



As anticipated above, we now fix the proportionality constant in ^'^'^'(r) ex g{r;lil2l) by dividing out the leading 
constants above so that the coefficient of 91i12o{t) is unity; e.g., g^^^{r) = .9220 (r) — 2.9221 {r) + 25222 {r). 
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X pljMC) pi{Kwc) r;(MC) r;(RHNC) zc(mc) 



1.0 


0.312 


0.317 


1.22 


1.23 


0.0526 


0.9 


0.309 


0.296 


1.05 


1.10 


0.0516 


0.8 


0.303 


0.277 


0.883 


0.949 


0.0487 


0.7 


0.287 


0.262 


0.714 


0.750 


0.0432 


0.6 


0.262 


0.254 


0.555 


0.562 


0.0350 


0.5 


0.234 




0.423 




0.0202 


0.4 


0.206 




0.333 




0.0202 


0.3 


0.175 




0.257 




0.0155 



TABLE I: Comparison of the estimated location of critical points from MC and RHNC data (see text). In the last column 
Zc = e^^ is the critical activity. 
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U/Ne 




/3/i 


I3P/P 


P (dP/dp)^ 


(to/ct 
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1.0 


-5.46 


-2.56 


-3.28 


0.64 


10.33 


1.031 


10.92 


0.9 


-4.88 


-1.81 


-1.79 


1.37 


11.26 


1.026 


9.76 


0.8 


-4.10 


-0.98 


-0.10 


2.24 


12.39 


1.020 


8.20 


0.7 


-3.25 


-0.19 


1.53 


3.08 


13.54 


1.014 


6.50 


0.6 


-2.49 


0.49 


2.88 


3.75 


14.48 


1.010 


5.00 


0.5 


-1.88 


1.05 


3.94 


4.25 


15.21 


1.007 


3.76 


0.4 


-1.26 


1.61 


5.00 


4.74 


15.95 


1.005 


2.51 


0.3 


-0.80 


2.00 


5.75 


5.10 


16.51 


1.003 


1.60 


0.2 


-0.39 


2.33 


6.39 


5.42 


17.00 


1.001 


0.79 


0.1 


-0.11 


2.54 


6.79 


5.62 


17.31 


1.000 


0.23 


0.0 


0.00 


2.61 


6.95 


5.69 


17.43 


1.000 


0.00 



TABLE IL Values of reduced internal and excess free energies, chemical potential, pressure, and inverse compressibility as a 
function of the coverage x for a fixed state point, p* — 0.7 and T* = 1.0. The last two columns report the reference HS diameter 
CTo (in units of a) and the average coordination number z. Expected errors are in the last digits. 
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FIG. 1: The two-patch Kern-Frenkel potential. Each sphere is divided into an attractive part (color code: green) and a repulsive 
part (color code: red). The attractive part is positioned on two symmetrically distributed patches identified by unit vectors 
n''' = fii and n'*'' = —hi {i — 1, 2), where the orientation vectors fii, n2 define angles 6i, 62 with the vector ri2 joining the 
centers of the two spheres and directed from sphere 1 to sphere 2. The particular case shown corresponds to a 40% fraction of 
attractive surface (coverage x)- 
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FIG. 2: Fluid-fluid coexistence lines of the two-patch model for different values of the coverage ranging from a full square-well 
potential down to x = 0.3. Points represent MC results while thick solid lines report RHNC values (for x ^ 0.6). Thinner solid 
lines are a guide for the eye whereas the thick dashed line shows the estimated GCMC critical point for any fixed coverage x- 
The MC data for x ~ 1-0 (SW) coincide within the numerical error with the ones reported in Refs. [s^ andfsil. 
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FIG. 3: Representation of the structures observed at small coverages x- (a-) The case of coverages such that each patch can 
be involved in only one interaction. In this case, the system forms polydisperse chains. The snapshot here refers to the case 
p* = 0.01 and T* = 0.07. (b) Values of x such that each patch can be involved in only two interactions. In this case, at low 
T the system forms bonded planes interacting with each other only via excluded volume interactions. The snapshot shows one 
such plane, (c) Values of x such that each patch can be involved in only three interactions. The crystal is now formed by 
interconnected planes, with a triangular arrangement of the particles in the plane. Adjacent planes are shifted in such a way 
that each particle sits in correspondence to the center of a triangle of the previous and following planes, (d) Values of x such 
that each patch can be involved in at most four interactions. The crystal is now formed by interconnected planes, with a square 
arrangement of the particles in the plane. Adjacent planes are shifted in such a way that each particle sits in correspondence 
to the center of a square of the previous and following planes. 
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FIG. 4: Behavior of ,g(12) from the RHNC equation for different coverages and two specific orientations of the patches: || 
configuration corresponding to fii • n2 ~ ±1 (a) and X configuration corresponding to fii • n2 = (b). All curves are for a 
state point with A = 1.5, p* = 0.7, and T* — 1.0. Black and green lines show the limiting cases of square-well (x = 1) and 
hard-sphere (x = 0) potentials, respectively. Other coverages are 0.9 — 0.1 for both the || and X configurations. 
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FIG. 5: Behavior of the averaged g{r, cos 6) for different coverages and two specific orientations of the patches: || configuration 
corresponding to fii • n2 = ±1 (a) and X configuration corresponding to ni • n2 = (b). All curves are for a state point with 
A = 1.5, p* = 0.7, and T* = 1.0. Both RHNC and MC results are depicted for x = 0, 0.1, 0.4, 0.5, 0.6, 0.9, 1.0. 





FIG. 6: Behavior of the structure factor Sooo{k) ((a) and (c) ) and the radial distribution function (7ooo(^) ((b) and (d)) for 
coverages x = 0.8 and x ~ 0-2, respectively. Here T* = 1.0 and A = 1.5 as before, while densities are p* = 0.7 and p* = 0.1. 
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FIG. 7: Angular distributions g{6, 62) as functions of cos d for two different orientations of the patches on sphere 2, given that 
sphere 1 is fixed with patches along the z axis. Results are reported for two different coverages, x = 0.8 ((a) and (b)) and 
X ~ 0.2 ((c) and (d)), and two different densities, p* = 0.7 ((a) and (c)) and p* = 0.1 ((b) and (d)), at the same temperature 
T* = 1.0. Again, the square-well width was set to A = 1.5. The colored arrows are cartoons of the orientations of sphere 2 
patches, corresponding to 62 — 0, 7r/2. Note that these are the same state points considered in Fig. (6] 
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FIG. 8: Plot of rotational invariants /i^^*^(r), /i''^^(r), and h^^^{r) as functions of r/a. Results are reported for two different 
coverages, x = 0.8 ((a) and (b) ) and x = 0.2 ((c) and (d)), and two different densities, p* = 0.7 ((a) and (c)) and p* = O.f ((b) 
and (d)), at the same temperature T* = 1.0 and square- well width A = 1.5. Again, these are the same state points considered 
in Figs. [n]and[7] 
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cos 6 cos 9 

FIG. 10: Angular distribution g{9, 62) for the one-patch model as a function of cos 9 for three different orientations of the patch 
on sphere 2, given that sphere 1 is fixed with patch along the z axis. This is the one-patch counterpart of Fig. [T] Results are 
reported for two different coverages, x = 0.8 ((a) and (b)) and x = 0.2 ((c) and (d)), and two different densities, p* = 0.7 ((a) 
and (c)) and p* = 0.1 ((b) and (d)), at the same temperature T* = 1.0 and square-well width A = 1.5. The colored arrows are 
cartoons of the orientations of the sphere 2 patch, corresponding to 62 ~ 0, 7r/2, n. 
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FIG. 11: Plot of rotational invariants h^^°{r), /i"^(r), and h^^°{r) for the one-patch model as functions of r/a. Results are 
reported for two different coverages, x = 0.8 ((a) and (b)) and x = 0-2 ((c) and (d)), and two different densities, p* = 0.7 ((a) 
and (c)) and p* = 0.1 ((b) and (d)), at the same temperature T* = 1.0 and square-well width A = 1.5. Again, these are the 
same state points considered in Fig. |5] 



